#ifndef _RECT_CIRCLEEXCLUDED_H_
#define _RECT_CIRCLEEXCLUDED_H_

#include "Core/Tensor.h"
#include "RectDomain.H"
#include "Circle.H"
template <int Dim>
class ExcludedCircle;

/**
 * The structure that describes the embedding of the irregular domain into the Cartesian grid.
 */
template <>
class ExcludedCircle<2> : public RectDomain<2>
{
public:
    enum
    {
        Dim = 2
    };
    using BaseClass = RectDomain<Dim>;
    template <class T>
    using vector = std::vector<T>;
    using rVec = Vec<Real, Dim>;
    using iVec = Vec<int, Dim>;

public:
    ExcludedCircle(const RectDomain<2> &aRectDomain,
                        const OrientedCircle &aCircle, Real atol);

public:
    const OrientedCircle embeddingCircle;
    const Real tol;
    Tensor<bool, Dim> inDomain;
    const Box<2> getMiniBox() const
    {
         return Box<2>{(iVec)floor(embeddingCircle.minimalEmbededdSquare()[0] / dx), 
         (iVec)floor(embeddingCircle.minimalEmbededdSquare()[1] / dx) + 1};
    }
    ///<这个数组的第i个代表unit(i)=const时的截断比矩阵。
    ///<第d个Tensor的(i,j)坐标存储了从（i,j）沿e_d正方向出发到intersectionPoint的相对网格距离比例.若0,则说明无交点。
    ///<这里首先假定每条mesh segment至多与圆只有一个交点。
    ///<于是我们可以把每个交点向下/左投影到最近的网格节点，并在这些网格存储对应线段的截断比\alpha.
    ///<其中\alpha代表了节点（关于下/左）的相对比例。
    Tensor<Real, Dim> cutRadios[2];
    const rVec cutPointMap(const iVec &idx, int d) const
    {
        return (idx + iVec::unit(d) * cutRadios[d](idx)) * dx;
    }
};

ExcludedCircle<2>::ExcludedCircle(const RectDomain<2> &aRectDomain,
                                            const OrientedCircle &aCircle, Real atol)
    : BaseClass(aRectDomain), embeddingCircle{aCircle}, tol(atol)
{
    vector<rVec> miniEmbededSquare = aCircle.minimalEmbededdSquare();
    const Box<2> miniBox = getMiniBox();
    Box<2> bx = *this;
    assert(bx.contain(miniBox));

    /// Step 0 Initialization.
    inDomain.resize(bx);
    inDomain = true;
    for (int i : {0, 1})
        cutRadios[i].resize(bx);

    /// Step 1 Determine nodeLabel and intersectionPoints.
    loop_box_2(miniBox, i, j)
    {
        iVec idx{i, j};
        rVec x = idx * dx;
        if (!aCircle.inEnclosure(x))
            inDomain(i, j) = false;
    }

    loop_box_2(miniBox, i, j)
    {
        iVec idx{i, j};
        rVec x = idx * dx;
        for (int d : {0, 1})
        {
            if (inDomain(i, j) != inDomain(idx + iVec::unit(d)))
            {
                const rVec CutP = aCircle.intersection2CartesianLine(1 - d, x);
                cutRadios[d](i, j) = (CutP[d] - x[d]) / (dx[d]);
            }
        }
    }
};
#endif // EMBEDEDGRID_H
